% rc9 Feb-2011 - 
%  data quality checks - b_norm, e_zero, esum, e_avg & e_03d
%
% $Id: e_o3d.m 492 2006-05-22 16:26:17Z edmaclin $
% e_o3d.m function version of eros_o3d.m - writes .avg files in opt-cont format
% - modified from optcont.m to use real SD info from eros_avg2.m 12/16/02
%
% 2004/05/10 declared SE arrays to avoid missing (n=0) bins
%            process only "good" channels
% 2004/05/19 individual .gch for each montage
% 2004/05/19 use .gch not .dis
% 2010/05/17 decimates if hdr.decimate > 1
% 2011/10/17 'good_chs' defined here to reflect "quality" summing...
% must have 'hdr.n_threshold' trials to be a 'good' channel

function e_o3d(in_path,out_path,log_file)

global hdr file_name i_block dc ac ph adc aac aph n2 ssdc ssac ssph ...
    msg zap_flg final_good_chs id_tags i_config min_n avg_PathName

% in_path = [hdr.data_path hdr.exp_name hdr.data_type hdr.hi '-' hdr.lo '\'];
% out_path = [hdr.data_path hdr.exp_name hdr.data_type hdr.hi '-' hdr.lo '\'];

disp('e_o3d')

e_o3d_id = '$Id: e_o3d.m 492 2006-05-22 16:26:17Z edmaclin $';
check = strcmp(e_o3d_id,id_tags); % Make sure we don't already have the tag
if sum(check) == 0
    id_tags{size(id_tags,1)+1,1} = e_o3d_id;  % Add it in then.
end

% Determine whether input and output files exist
%  !!! 'exist' function doesn't work in compiler :-(
%
k=findstr('.',file_name);
in_file_name=[file_name(1:k(1)) 'avm'];
in_file_flg=fopen([avg_PathName '\' in_file_name]);
out_file_name=[file_name(1:k(1)) 'avg'];
out_file_flg=fopen([avg_PathName '\' out_file_name]);
fclose('all');

% Open the log file if it exists
%
if ~strcmp(log_file,'')
    logid = fopen(log_file,'a');
else logid = 1; end

% fprintf(logid,'%s flag=%2.0f\n',[in_path in_file_name],in_file_flg)
% fprintf(logid,'%s flag=%2.0f\n',[out_path out_file_name],out_file_flg)

% if (exist([in_path in_file_name]) && ~exist([out_path out_file_name]))
if in_file_flg > 0 && out_file_flg < 0
    
    % see if hdr.decimate exists & set to 1 if not
    %
    try
        tmp=hdr.decimate;   % 78 = 2/sec for ECC @39.0625 Hz
    catch
        hdr.decimate=1;    % apparently don't want to decimate
    end
    
    fprintf(logid,'Loading %s\n',[avg_PathName '\' in_file_name]);
%     load([in_path in_file_name],'adc','aac','aph','n2','SDdc','SDac','SDph','mrdc','mrac','-mat'); 
    load([avg_PathName '\' in_file_name],'-mat');
    disp('Loaded')
%     disp(['size(SDac) = ' num2str(size(SDac))])
    
    % decimate (all channels) if desired
    %
    if hdr.decimate > 1
        disp(['Decimating by ' num2str(hdr.decimate)])
        for i_ch=1:size(adc,2)
            for i_cond=1:size(adc,3)
%                 disp(['i_ch = ' num2str(i_ch) ';i_cond= ' num2str(i_cond)])
                aac_tmp(:,i_ch,i_cond)=decimate(aac(:,i_ch,i_cond),hdr.decimate);
                adc_tmp(:,i_ch,i_cond)=decimate(adc(:,i_ch,i_cond),hdr.decimate);
                aph_tmp(:,i_ch,i_cond)=decimate(aph(:,i_ch,i_cond),hdr.decimate);
                SDac_tmp(:,i_ch,i_cond)=decimate(SDac(:,i_ch,i_cond),hdr.decimate);
                SDdc_tmp(:,i_ch,i_cond)=decimate(SDdc(:,i_ch,i_cond),hdr.decimate);
                SDph_tmp(:,i_ch,i_cond)=decimate(SDph(:,i_ch,i_cond),hdr.decimate);
            end
        end
%         disp([num2str(size(aac_tmp,1)) ' points'])
        
        clear aac
        aac=aac_tmp;
        clear adc
        adc=adc_tmp;
        clear aph
        aph=aph_tmp;
        clear SDac
        SDac=SDac_tmp;
        clear SDdc
        SDdc=SDdc_tmp;
        clear SDph
        SDph=SDph_tmp;
    else
         disp('not decimating')        
    end % decimate
   hdr.decimate=1;  % !!!
    time=(1:size(adc,1))-fix(hdr.avg_baseline/hdr.decimate);
    time=time*hdr.period_ms*hdr.decimate;
    
    % adjust time to seconds for pretty display in opt3d (doesn't work?)
    %
%     if time(end)-time(1) > 1
%         time=time/1000;
%     end
    
    adc=adc+1.0;     % opt3d needs mean = 1
    aac=aac+1.0; 
    
    fprintf(logid,'Opening file %s\n',[avg_PathName '\' out_file_name]);
    fid=fopen([avg_PathName '\' out_file_name],'w');
    
    fprintf(logid,'Converting...\n');    
  
    % define good channels as those with > n_thresh trials in EVERY bin
    %
    min_n=min(n2,[],2); % n of smallest bin in each channel
    disp(['attempting to define final_good_chs; size(f_g_c) = ' num2str(size(final_good_chs))])
    disp(['i_config = ' num2str(i_config)])
    final_good_chs{i_config}=find(min_n>=hdr.n_threshold);
%     good_chs(i_config)=find(hdr.n_good.tot>=hdr.n_threshold);
    
    % print mean values to .avg file
    %
    for j_chan=1:length(final_good_chs{i_config})
        i_chan=final_good_chs{i_config}(j_chan);
        fprintf(fid,'CHANNEL %4.0f,',i_chan);
        fprintf(fid,'%12.0f,',n2(i_chan,:),n2(i_chan,:),n2(i_chan,:));
        fprintf(fid,'\r\n');
        fprintf(fid,'           0,');
        fprintf(fid,'%12.0f,',n2(i_chan,:),n2(i_chan,:),n2(i_chan,:));
        fprintf(fid,'\r\n');
        for i_time =1:size(adc,1)
            
            fprintf(fid,'%12.5f,',time(i_time));
            
            for i_cond=1:hdr.n_conds
                fprintf(fid,'%12.5f,',adc(i_time,i_chan,i_cond));
            end
            
            for i_cond=1:hdr.n_conds
                fprintf(fid,'%12.5f,',aph(i_time,i_chan,i_cond));
            end
            
            for i_cond=1:hdr.n_conds                
                fprintf(fid,'%12.5f,',aac(i_time,i_chan,i_cond));
            end

            fprintf(fid,'\r\n');
            
        end % i_time
    end % i_chan

    SEdc=zeros(size(SDdc));
    SEac=zeros(size(SDac));
    SEph=zeros(size(SDph));
%     disp(['size(SEdc) = ' num2str(size(SEdc))])

    % Convert Standard Deviation to Std Error
    %
    for i_cond=1:hdr.n_conds
        for j_chan=1:length(final_good_chs{i_config})
            i_chan=final_good_chs{i_config}(j_chan);
%             fprintf('i_chan=%3.0f, n2=%3.0f, sqrt(n2)=%3.0f\n',i_chan,n2(i_chan,i_cond),sqrt(n2(i_chan,i_cond)))

            if n2(i_chan,i_cond) > 0
                SEdc(:,i_chan,i_cond)=SDdc(:,i_chan,i_cond)/sqrt(n2(i_chan,i_cond));
                SEac(:,i_chan,i_cond)=SDac(:,i_chan,i_cond)/sqrt(n2(i_chan,i_cond));
                SEph(:,i_chan,i_cond)=SDph(:,i_chan,i_cond)/sqrt(n2(i_chan,i_cond));
            else
                SEdc(:,i_chan,i_cond)=0.0001;
                SEac(:,i_chan,i_cond)=0.0001;
                SEph(:,i_chan,i_cond)=0.0001;
            end
        end % i_chan
    end % i_cond
%     SEdc=SDdc;SEac=SDac;SEph=SDph;
    
    % Print the Std Error
    %
    for j_chan=1:length(final_good_chs{i_config})
        i_chan=final_good_chs{i_config}(j_chan);
        fprintf(fid,'CHANNEL %4.0f,',i_chan);
        fprintf(fid,'%12.0f,',n2(i_chan,:),n2(i_chan,:),n2(i_chan,:));
        fprintf(fid,'\r\n');
        fprintf(fid,'           0,');
        fprintf(fid,'%12.0f,',n2(i_chan,:),n2(i_chan,:),n2(i_chan,:));
        fprintf(fid,'\r\n');
        for i_time =1:size(adc,1)
            
            fprintf(fid,'%12.5f,',time(i_time));

                           for i_cond=1:hdr.n_conds
                                if SEdc(i_time,i_chan,i_cond) >= 0.0001
                                    fprintf(fid,'%12.5g,',SEdc(i_time,i_chan,i_cond));
                                else
                                    fprintf(fid,'%12.5g,',0.0001);
                                end
                            end

                           for i_cond=1:hdr.n_conds
                                if SEph(i_time,i_chan,i_cond) >= 0.0001
                                    fprintf(fid,'%12.5g,',SEph(i_time,i_chan,i_cond));
                                else
                                    fprintf(fid,'%12.5g,',0.0001);
                                end
                            end

                           for i_cond=1:hdr.n_conds
                                if SEac(i_time,i_chan,i_cond) >= 0.0001
                                    fprintf(fid,'%12.5g,',SEac(i_time,i_chan,i_cond));
                                else
                                    fprintf(fid,'%12.5g,',0.0001);
                                end
                            end
                           
            fprintf(fid,'\r\n');
            
        end % i_time
    end % i_chan
    msg=sprintf('Writing %s\n',[avg_PathName '\' out_file_name]);
    fprintf(logid,'Writing %s\n',[avg_PathName '\' out_file_name]);
    fclose(fid); 
else % if various files exist or don't...
    if in_file_flg < 0
        msg=sprintf('%s does not exist.\n',[avg_PathName '\' in_file_name]);
        fprintf(logid,'%s does not exist.\n',[avg_PathName '\' in_file_name]);
    end
    if out_file_flg > 0
        msg=sprintf('%s already exists!\n',[avg_PathName '\' out_file_name]);
        fprintf(logid,'%s already exists!\n',[avg_PathName '\' out_file_name]);
    end
end

if logid ~= 1, fclose(logid); end